Chapter 2 



Integrals as sums and 
derivatives as differences 



We now switch to the simplest methods for integrating or differentiating a 
function from its function samples. A careful study of Taylor expansions 
reveals how accurate the constructions are. 

2.1 Numerical integration 

Consider a function / - we'll specify which assumptions we need to make 
about it in a minute. Let us reformulate the integral 



by breaking up the interval [a, b] into subintervals [ ], with Xj = jh, 

h=l/N, and 

= Xo < X\ < . . . < Xn — 1- 

Together, the (xj) are called a Cartesian grid. Still without approximation, 
we can write 



A quadrature is obtained when integrals are approximated by quantities that 
depend only on the samples f(xj) of / on the grid Xj. For instance, we can 
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approximate the integral over [xj-i, Xj] by the signed area of the rectangle 
of height f(xj-i) and width h: 




f(x)dx ~ hf(x j ^ 1 ). 



Putting the terms back together, we obtain the rectangle method: 

„1 N 

/ f(x)dx ~ h V/(^_i). 
Jo j=1 

(Insert picture here) 

To understand the accuracy of this approximation, we need a detour 
through Taylor expansions. Given a very smooth function f(x), it sometimes 
makes sense to write an expansion around the point x = a, as 

/(*) = f(a) + f'(a)(x -a) + \f"(a)(x - af + ^f"(a)(x - a) 3 + . . . 
As a compact series, this would be written as 

oo 1 

f{x) = Y.-/ n) ^ x - a ) n ' 

n=0 

where n\ — n • (n — 1) • . . . 2 • 1 is the factorial of n. The word "sometimes" is 
important: the formula above requires the function to be infinitely differen- 
tiable, and even still, we would need to make sure that the series converges 
(by means of the analysis of the previous chapter). The radius of conver- 
gence of the Taylor expansion is the largest R such that, for all x such that 
\x — a\ < R, the series converges. 

It is also possible to truncate a Taylor series after N terms, as 

N-l 

= E n\ fin){a){x ~ ar + m fiN) ^ x - a)7V ' (2 - 1} 

n=0 

where y is some number between x and a. The formula does not specify what 
the number y is; only that there exists one such that the remainder takes 
the form of the last term. In spite of the fact that y is unspecified, this is a 
much more useful form of Taylor expansion than the infinite series: 
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• We often only care about the size of the remainder, and write inequal- 
ities such as 

1 f (N \y)(x-a) N \<±- ( max\fW(z)\) \x - a\ N , 



N\ y /v N\ \ze[x,a] 

or, for short, 

^ N \y)(x-af = 0(\x-an 

where all the "constants" are hidden in the 0. So, the remainder is 
on the order of \x — a\ N . The limit implicit in the notation is here 
x — a — > 0. 

• The finite- N formula is valid whenever the function is N times differ- 
ent iable, not infinitely differentiable! In fact, the Taylor series itself 
may diverge, but equation (2.1) is still valid. 

In order to study the problem of approximating J^ 3 i f(x) dx, let us ex- 
pand / in a (short, truncated) Taylor series near x = Xj-i: 

f(x) = f(x^) + f'(y(x))h, 

where y(x) G [ ]. We've written y(x) to highlight that it depends on 
x. This works as long as / has one derivative. Integrating on both sides, 
and recognizing that f(xj-i) is constant with respect to x, we get (recall 

x j x j_i — hi j 

nXj PXj 

/ f(x) dx = hf(xj-i) +h f'(y(x)) dx. 

J Xj — l J Xj — i 

We don't know much about y(x), but we can certainly write the inequality 

nXj PXj 

h \ / f'(y(x)) dx\ < h max \f'(y)\ / ldx = h 2 max \f'(y)\. 
Jxj-! ye[xj-i,xj] Jx 3 -i v£[xj-i,xj] 

So, as long as / is differentiable, we have 

f(x) dx = hf(x^ 1 )+0(h 2 ), 



r 

Jxj-1 
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where the derivative of / hides in the sign. The integral over [0, 1] has iV 
such terms, so when they are added up the error compounds to N • 0(h 2 ) = 
0(h) (because h = 1/-/V). Thus we have just proved that 

.1 N 

/ f(x)dx~hJ2f(xj-i) + 0(h). 

The error is on the order of the grid spacing h itself, so we say the method 
is first-order accurate (because the exponent of h is one.) 

Choosing the (signed) height of each rectangle from the left endpoint Xj-i 
does not sound like the most accurate thing to do. Evaluating the function 
instead at £j-i/2 = (j — 1/2) /i is a better idea, called the midpoint method. 
It is possible to show that 

/ f(x)dx^hJ2f(xj-i/2) + 0(h 2 ), 
Jo j=1 

provided the function / is twice different iable (because /" hides in the O 
sign). The accuracy is improved, since h 2 gets much smaller than h as h — > 0. 
We say the midpoint method is second-order accurate. 

Another solution to getting order-2 accuracy is to consider trapezoids 
instead of rectangles for the interval [ ]. The area of the trapezoid 

spanned by the 4 points 0); (xj-i, f(xj-i)); (xj, f(xj)); (xj, 0) is h(f(xj-i)+ 

f(xj))/2. This gives rise to the so-called trapezoidal method, or trapezoidal 
rule. 

(Insert picture here) 

Let us compare this quantity to the integral of / over the same interval. 
Consider the truncated Taylor expansions 

f{ X ) = /(xj-!) + ffa-Mx - Xj-i) + 0(h 2 ), 

f(x,) = f(x,. 1 ) + f'(x^ 1 )h + 0(h 2 ), 

where the second derivative of / appears in the O sign. Integrating the first 
relation gives 

P f(x) dx = hf(x^) + ^f'(x^) + 0(h 3 ). 

Jxj-l A 
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The area of the trapezoid, on the other hand, is (using the second Taylor 
relation) 

+ /(*;)) = hf{x^) + y/'fo-i) + 0(/* 3 ). 

Those two equations agree, except for the terms 0(h 3 ) which usually differ 
in the two expressions. Hence 



Jxi-l 



f( x )dx = ^(f(x J _ 1 ) + f(x J )) + 0(h 3 ). 



We can now sum over j (notice that all the terms appear twice, except 
for the endpoints) to obtain 



N-l 



f(x) dx = |/(0) + h J2 /(*;) + + 0(h 2 ). 



3=1 



The trapezoidal rule is second-order accurate. All it took is a modification 
of the end terms to obtain 0(h 2 ) accuracy in place of 0(h). 
Example: f(x) = x 2 in [0,1]. We have Q \ x 2 dx = 1/3. 

• For the rectangle rule with N = 4 and h — 1/4, consider the gridpoints 
x = 0,1/4,1/2,3/4, and 1. Each rectangle has height f(xj-i) where 
Xj-i is the left endpoint. We have 



1 2, 1 

x dx ~ - 







1 1 3 21 

+ 4 2 + 2 2 + 4 2 



14 

64 



.2188... 



This is quite far (0(h), as we know) from 1/3. 

For the trapezoidal rule, consider the same grid. We now also consider 
the grid point at x — 1, but the contribution of both x = and x — 1 
is halved. 



x dx 



11 1 3 2 1 ' 
2- + 4^ + 2^ + 4^ + 2- 1 



22 
64 



.3438... 



This is much closer (0(h 2 ) as a matter of fact) from 1/3. 

We'll return to the topic of numerical integration later, after we cover 
polynomial interpolation. 
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2.2 Numerical differentiation 

The simplest idea for computing an approximation to the derivative u'(xj) 
of a function u from the samples Uj = u(xj) is to form finite difference ratios. 
On an equispaced grid Xj = jh, the natural candidates are: 

• the one-sided, forward and backward differences 

U j+1 - Uj _ Uj - Uj-x 

A+Uj - h ' A - " h ' 

• and the two-sided centered difference 

U j+1 - 



A Uj = 



2h 



Let us justify the accuracy of these difference quotients at approximating 
derivatives, as h — > 0. The analysis will depend on the smoothness of the 
underlying function u. 

Assume without loss of generality that Xj-i = —h, Xj = 0, and Xj+i = h. 
For the forward difference, we can use a Taylor expansion about zero to get 

u (h) = u(0) + W(0) + e e [0, h\. 

(The greek letter ^ is "xi".) When substituted in the formula for the forward 
difference, we get 

u(h)-u(0) =AQ) + h u „ { ^ 

The error is a 0(h) as soon as the function has two bounded derivatives. We 
say that the forward difference is a method of order one. The analysis for 
the backward difference is very similar. 

For the centered difference, we now use two Taylor expansions about zero, 
for u(h) and u(—h), that we write up to order 3: 

u(±h) = u(0) ± hu'(0) + ^«"(0) ± 

with ^ either in [0, h] or [—h, 0] depending on the choice of sign. Subtract 
u(—h) from u(h) to get (the h 2 terms cancel out) 
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The error is now 0(h 2 ) provided u is three times differentiable, hence the 
method is of order 2. 

The simplest choice for the second derivative is the centered second dif- 
ference 

2 _ Uj+l - 2Uj + Uj-i 

It turns out that A 2 = A + A_ = A_A + , i.e., the three-point formula for the 
second difference can be obtained by calculating the forward difference of the 
backward difference, or vice-versa. A 2 is second-order accurate: the error is 
0(h 2 ). To see this, write the Taylor expansions around x = to fourth order: 

u(±h) = «(0) ± hu'(0) + ^«"(0) ± ju"'(0) + ^""(0- 

with £ either in [0, h] or [—h, 0] depending on the choice of sign. The odd 
terms all cancel out when calculating the second difference, and what is left 
is 

^•+i- 2 ^ + ^--i = u // (0 ) + ^ u ""(0- 

ft \-Zi 

So the method really is 0(h 2 ) only when u is four times differentiable. 
Example: f(x) = x 2 again. We get f'(x) = 2x and f"(x) = 2. 

. . (x + h) 2 -x 2 2xh + h 2 
A+/(x) = '- = = 2x + h. 

A_/(x) = 



(x + h) 2 


- x 2 2xh + h 2 


h 


h 


x 2 — (x - 


- h) 2 2xh - h 2 


h 


h 


_ (x + h) 


2 -(x- h) 2 Axh 




2h 2h 


(x + h) 2 - 


- 2x 2 + (x- h) 2 



A of(x) = ± '—^ '- = — = 2x. 



2^ 2 

h 2 h 2 
The error is manifestly h for the forward difference, and —h for the backward 
difference (both are 0(h).) Coincidentally, the error is zero for A and A 2 : 
centered differences differentiate parabolas exactly. This is an exception: 
finite differences are never exact in general. Of course, = 0(h 2 ) so there's 
no contradiction. 

When u has less differentiability than is required by the remainder in the 
Taylor expansions, it may happen that the difference schemes may have lower 
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order than advertised. This happens in numerical simulations of some differ- 
ential equations where the solution is discontinuous, e.g., in fluid dynamics 
(interfaces between phases). 

Sometimes, points on the left or on the right of Xj are not available 
because we are at one of the edges of the grid. In that case one should use 
a one-sided difference formula (such as the backward or forward difference.) 
You can find tables of high-order one-sided formulas for the first, second, and 
higher derivatives online. 

The more systematic way of deriving finite difference formulas is by dif- 
ferentiating a polynomial interpolant. We return to the topics of integration 
and differentiation in the next chapter. 
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